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c+h ' Abstract 

O' 

We introduce a model of granular matter and use a stress ensemble to analyze shearing. 

Monte Carlo simulation shows the model to exhibit a second order phase transition, 

associated with the onset of dilatancy. 
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1. Introduction. 

Static granular matter, such as a sand pile, can exist in a range of densities. In its 
densely packed state, the more common form, it expands under shear, while when its 
grains are loosely packed it contracts under shear [RN]. We introduce a model within 
which the transition between these qualitative behaviors is singular in the precise 
traditional sense of a second order phase transition. (For two dimensional treatments 
see [AS, PLR], and for an interpretation as a glass transition see [CH].) 

Our "granular hard cubes" model is a granular variant of the classical hard cubes 
model of equilibrium statistical mechanics [HHB] , the latter being a simplification of 
the classical hard sphere model in which the spheres are replaced by nonoverlapping, 
parallel unit cubes. In our granular version, following Edwards [EO], we only allow 
configurations in which the parallel cubes are mechanically stable under gravity, where 
gravity is acting parallel to edges of the cubes. We require that the centers of our cubes 
be approximately at the vertices of an fee lattice, which we orient so that gravity acts 
along the 001 axis, so the fee lattice is thought of as parallel 2-dimensional square 
lattice layers. We enforce this approximate structure by requiring that each cube sit 
on exactly four cubes in the layer below it; see Fig. 1. 




Figure 1. Layers of parallel cubes. 

We want to model the reaction of a sandpile to shear and in particular we need 
a model which exhibits a critical state density, a density separating a high density 
regime of states, in which the material expands under shear, from a low density regime 
in which it contracts under shear [RN] . As we show, this is easy to arrange. The more 
interesting question is whether the model predicts that the transition through that 
critical state density is smooth, or is singular. 
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To understand why it is natural for a granular model to have a critical state 
density it is useful to use the following "stress" ensemble. Instead of constant density 
we introduce a pressure parameter, and in place of constant strain we introduce a 
shear stress parameter. In other words the states of the model, which, in the general 
approach of Edwards [EO] are probability densities of configurations of grains, optimize 
the free energy 

F{pJ):=S-pV + faV, (1) 

where S = ln(Z p j) is the entropy, V is the volume in physical space, and a is the 
shear strain. The parameter p is pressure, / is shear stress, and 



j pJ 




exp(— pV + faV) dVda. 



(2) 



This is equivalent to saying that a configuration C has an unnormalized probability 
density: 

Pr(C) = exp(-pV + faV). (3) 

It follows from this probability density that, in the absence of shear, the state of 
lowest possible particle density, called random loose packing, occurs at zero pressure 
[AR1] and some shear strain a = olq. Therefore at sufficiently low pressure if one 
shears the system from oto, the density has no way to change except to increase. So at 
sufficiently low density there is a regime in which the system contracts. 

On the other hand, at least since Reynolds [Re] one understands the expansion of 
(dense) granular matter under shear, which he termed dilatancy, in geometric terms, 
as the need for parallel layers of spheres to get out of each others' way under a strain 
deforming the layers, either by separating and/or by thinning within the layers. So we 
adjust the model to arrange for this (thinning) phenomenon, as follows. 



a) 




b) 




Figure 2. a) A view from above of a layer in a configuration, with a definition of the 
angle a; b) A cubic grain, with bumps on top and bottom faces. 

To react to the strain (see Fig. 2a) we mimic the spherical caps that grains in 
one layer present to grains in neighboring layers by adding "bumps" to the top and 



bottom of our cubes; see Fig. 2b. Such cubes will, when we strain a sufficiently dense 
configuration in response to the shear, feel the bumps in the layers above and below 
it, and this will produce a dilatancy effect by thinning the layers. By considering 
optimally symmetric configurations one quickly finds 4>d '■= (1 + w ) _1 5 where w is the 
width of a bump, as an estimate of the density required for dilatancy; see Fig. 3. In 
our simulations we use w = 0.15, and so 4>d ~ 0.59. The qualitative dilatancy effect is 
not sensitive to the precise value of w, though the value cf)d is- 

So the bumps, together with the universal behavior at low pressure or density, 
explain the existence of a critical state density. The main question then is: as density 
is varied, does the transition through the critical state density proceed in a smooth 
or in a singular manner? We will show there is an unambiguous second order phase 
transition at dilatancy onset, at density near (p^. This suggests, by analogy with 
matter in thermal equilibrium, that the material in the two regimes differs in other 
characteristics as well, for instance it would be expected that the yield force behave 
differently in the two regimes. In [SNRS] the yield force is measured as a function of 
density, and a second order phase transition is found at a density roughly 0.598. The 
critical state density was not carefully measured because of the experimental setup, but 
approximately coincided. Our result suggests the two phenomena are simply different 
manifestations of the same transition. 

We also note that the behavior of granular matter at the random close packing 
density, about 0.64, has been interpreted in [Ra, AR2] as a first order phase transition 
in which the high density regime is an ordered phase. Together with our results here, 
this means that the usual freezing transition of equilibrium fluids, at which the material 
acquires the solid features of an ordered (typically crystalline) internal structure as well 
as a strong resistance to shear, seems to be split into two stages in granular materials, 
with the resistance to shear occuring at dilatancy onset at density about 0.6 and the 
ordered internal structure occuring at the random close packing density of about 0.64. 
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Figure 3. A regular arrangement of grains at maximum density, from above, 
configuration can be partitioned into unit cells defined by the dotted lines. 
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2. Simulating the Model. 

As noted above, our model consists of hard, parallel, oriented unit cubic grains, 
with two small cubic bumps of width w < 1 attached to the middle of their upper 
and lower faces. We use w = 0.15. Each grain must sit on exactly four other grains 
(except grains on the outside of the configuration), but not on the bumps, and grains 
cannot overlap. Thus, the grains appear on discrete vertical levels, so the distance in 
the vertical or ^-direction between the centers of grains on adjacent levels is equal to 
1; see Fig. 4. 




Figure 4. An arrangement of grains, viewed from above. 

Each simulation begins with a perfectly regular arrangement of N = n 3 grains in 
a large cubic structure, so that there are n 2 grains in each level, with n in each row 
and each column of a level, and n total levels. (For ease in notation we associate a 
row with the x-direction and a column with the y-direction; note that the model has 
a well-defined row and column structure.) The grains on the side boundaries of the 
cubic structure are allowed to make coordinated moves corresponding to changes in 
configuration volume and shape, but they are not allowed to move singly. 

We allow three types of Monte Carlo steps. With probability 1 — 1/N we choose 
a random grain inside the boundary to move in either the x- or ^-direction; and 
with probability 1/N we allow a configuration to change shape or volume. When a 
configuration changes volume, the x- and y-coordinates of each grain, on each level, 
are scaled by the same factor A. When it changes shape, the ^-coordinate of each grain 



in the kth column of each level is scaled, for some a, by the factor ak. 

To prevent grains from getting stuck between the boundary grains, we impose 
top-bottom periodic boundary conditions, and we remove the bumps from the grains 
nearest to the boundary and on the boundary. 

As is standard in Monte Carlo, to simulate the distribution (3) we choose steps 
with relative acceptance probabilities given by V N exp(— pV + faV). (The factor V N 
arises from the necessity of rescaling position coordinates to live in a fixed-size box 
[Mc]). Here a is the angle of deformation of a configuration (see Fig. 2b), where 
from symmetry we can assume a > without loss of generality. The volume V of a 
configuration is the volume of the convex hull of the set of centers of all the grains. 
The volume can be computed directly from the scale factor A described above. 

3. Results of the Simulations. 

We ran Markov chain Monte Carlo simulations on the model, with the Monte 
Carlo steps designed so that the stationary distribution of the Markov chain has the 
probability density 

exp(— pV + faV) 



in , 



AC) 



j pj 



(4) 



described above. We want to measure how the (average) density <fi = N/V changes 
as shear stress / varies near / = 0. (For convenience we do not include the bumps in 
the calculation of </>, but since particle number N is fixed in our simulation, the "true" 
volume fraction is just a constant multiple of <p.) That is, we want to estimate the 
derivative 

d{</>) 



D 



df 



(5) 



/=o 



of the average of <p as a function of p. From the definition of m p t it is easy to see that 



D = D(p) = N(a) 



p.O 



(4>)p,o {Va) pfi 



(6) 



where {-) p j denotes average with respect to m p j. Thus, to estimate D(p) from our 
simulations we set / = and calculate the sample averages of Na, (f), and Va. With 
the free energy F(p, f) = S — pV + faV, which is the quantity minimized by the state 
Tn p j, one can see from a direct calculation that 



d(p) = -(Koii ( d -i 



df \ dp 



(7) 



/=o 



which shows that the function D(p) is a quantity for which a discontinuity might 
reasonably be described as a second order phase transition. 




Figure 5. The graph of D(p) 



for systems of 14 3 = 2744, 12 3 = 1728, 



/=o 



10 3 = 100, and 8 3 = 512 grains (datal-data4, resp.). The insert gives error bars on 
the system with 14 3 grains. 



We investigate systems of 8 3 = 512, 10 3 = 100, 12 3 = 1728 and 14 3 = 2744 grains 
at pressures p\ < ... < p^, where the pressures chosen are different for each size system. 
In the simulations we begin at low pressure p\, and then slowly increase p until just 
after the derivative D(p) falls sharply below zero, roughly corresponding to pressure 
Pk- The last configuration in the simulation of pi is used as the starting configuration 
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of the simulation of Pi+i- We use the standard biased autocorrelation function to 
determine a "mixing time" , measured as the number of Monte Carlo steps required 
before the autocorrelation first crosses zero. We run our simulations long enough so 
that the simulation of each pi contains, on average, at least 20 mixing times, and we 
run 200 independent copies of each simulation to obtain the averages (Na) Pj o, (4>)p,o, 
and (Va) Pj o. From these averages we compute D(p) from (6). (We also ran some of 
our simulations much longer, with fewer copies and pressures, and noted agreement 
with the already-obtained data on D(p).) Then, for error bars on D(p), we repeat the 
entire experiment 8 times and use the Student's t-distribution. 

Our simulations suggest that D(p) develops a discontinuity as the system size 
increases (see Fig. 5). Changing variables we note that D, as a function of (0) Pj o> 
develops a discontinuity near the (average) volume fraction 4>d pa 0.59 discussed in the 
introduction; we plot (0) p ,o against p in Fig. 6b and D against (4>)p,o in Fig. 6a. 

Note that for p < 1, D(p) exhibits regular behavior in which D(p) is roughly 
constant, D(p) = rj pa 0.001. Then for 1 < p < 2, D(p) has some oscillation which is 
characteristic to the system size. Then for p > 2, D(p) steadily increases to a peak, 
then sharply decreases through zero. We believe that D(p) is discontinuous in the 
limit of infinite system size, but the rate of change of D(p) is so large it is difficult 
to measure its variation with system size. Instead we consider two measures of the 
interval R in which the discontinuity is developing, and then note that R gets smaller 
and smaller as system size increases. 

We expect that the oscillation observed in D(p) in the interval 1 < p < 2 is 
caused by finite-size effects, so that it disappears in the limit N — y oo. Furthermore 
we expect the critical pressure p c to fall at either the end or the beginning of the 
oscillation region. These two possibilities underlie, respectively, our measurements R\ 
and i?2, defined below. Let r\ be the average value of D(p) over p < 1. The left 
endpoint of R\ is the smallest value of p where p > 2 and D(p) pa rj, and the right 
endpoint of R\ is the value of p where p > 2 and D(p) pa 0. The widths of R\ are 
approximately 1.75 ± 0.05, 1.15 ± 0.05, 0.75 ± 0.05, and 0.50 ± 0.05 for systems of 
N = 8 3 , 10 3 , 12 3 and 14 3 grains, respectively. On the other hand, the left endpoint 
of i?2 is defined as the smallest value of p such that D(p) differs significantly from r/, 
while the right endpoint of R2 is the same as the right endpoint of R\. We estimate 
that the widths of R2 are approximately 3.55 ± 0.05, 1.55 ± 0.05, 1.25 ± 0.15, and 
1.15 ± 0.05 for systems of N = 8 3 , 10 3 , 12 3 and 14 3 grains, respectively. 

In either case, based on the decreasing sizes of the intervals R\ and R2, we believe 
that the jump from D(p) pa v to D(p) <C occurs over an interval which is vanishing in 
the infinite volume limit, which, given the relation between D(p) and the free energy 
F(p, f) in (7), is reasonably termed a second order phase transition. And as noted in 
the introduction, there is experimental evidence for this interpretation in [SNRS]. 
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Figure 6. a) D(p) as a function of (<f>) p ,o', b) (0) p ,o as a function of p. In both plots, 
datal-data4 represent systems of 14 3 = 2744, 12 3 = 1728, 10 3 = 100, and 8 3 = 512 
grains, respectively. 
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